Quench dynamics and non equilibrium phase diagram of the Bose-Hubbard model 
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We investigate the time evolution of correlations in the Bose-Hubbard model following a quench 
from the superfluid to the Mott insulator. For large values of the final interaction strength the 
system approaches a distinctly non-equilibrium steady state that bears strong memory of the initial 
conditions. In contrast, when the final interaction strength is comparable to the hopping, the 
correlations are rather well approximated by those at thermal equilibrium. The existence of two 
distinct non-equilibrium regimes is surprising given the non-integrability of the Bose-Hubbard model. 
We relate this phenomenon to the role of quasi-particle interactions in the Mott insulator. 
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Recent experiments with ultra cold atomic gases have 
opened exciting possibilities for studying non-equilibrium 
quantum dynamics of many body systems. In particular, 
the high degree of tunability allows one to rapidly change 
system parameters and observe the subsequent quantum 
evolution. Furthermore, thanks to the almost perfect iso- 
lation of the atoms from the environment, the quantum 
dynamics can remain coherent for exceedingly long times. 
These advantages were used, for example, to study non 
adiabatic dynamics across the quantum phase transition 
between a superfluid and a Mott insulator Q, Q , as well 
as the crossover of paired fermion superfluids from weak 
to strong coupling [H, 0] . 

In many cases the system parameters are changed so 
fast, that one may consider the sudden limit: The system 
is prepared in the ground state of an initial Hamiltonian 
Hi, and then evolves under the influence of a different 
Hamiltonian Ht. Fundamental questions that arise con- 
cern the approach of the system to a new steady state 
and the nature of this steady state. Does it retain mem- 
ory of the initial state? How is it related to the thermal 
equilibrium of Hp. These questions were addressed in 
a number of recent works, by solving various integrable 
models [H 0, 0, @, H El ■ In these systems the long time 
steady state was found to be non thermal and often car- 
ried memory of the initial state. In a fascinating exper- 
iment, Kinoshita et al. investigated the thermaliza- 
tion of strongly interacting ultra cold atoms in a nearly 
integrable situation. The result, at the maximal time 
scale of the experiment, was a non thermalized steady 
state. 

In this Letter we numerically investigate the evolution 
of correlations following a sudden change of parameters 
in the non-integrable Bose Hubbard model fBHM) |l 2ll. 
describing cold atomic gases in optical lattices [l3|, Il4| 
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In equilibrium at integer filling, this model exhibits a 



quantum phase transition at a critical value of the inter- 
action strength U/J= u c , between a superfluid (U/J < 
u c ) and a Mott insulating state (U/J > u c ) [12j. In a 
one-dimensional system with unit filling v} c D ~ 3.37 (ni | 
and in two dimensions u 2 c D w 16.7 (l6| . 

Our study is motivated by the experiment with ultra 
cold bosonic atoms in an optical lattice 0] , where the lat- 
tice intensity was increased suddenly, taking the system 
from the superfluid phase into the Mott insulator regime. 
Following the quench, a remarkable series of collapse and 
revivals of the interference pattern was observed, which 
relaxed after a few oscillations. What processes are re- 
sponsible for the relaxation and what is the nature of 
the steady state that is reached? The general expecta- 
tion for non-integrable models like |T]) is that the long 
time steady state will be essentially equivalent to ther- 
mal equilibrium. Surprisingly, we find that this is not 
always the case. When the interaction Uf in the final 
state is much larger than J the system reaches a quasi- 
steady state that is very different from thermal equilib- 
rium and retains memory of the initial state. As Uf de- 
creases the nature of the steady state changes, and in 
the region where Uf is comparable to J the steady state 
correlations are well approximated by those at thermal 
equilibrium (cf. Fig. [4]). We relate this crossover in the 
non-equilibrium behavior to the ineffectiveness of quasi- 
particle interactions in the strongly interacting regime. 

The revivals seen in the interference patterns in the ex- 
periment are easy to understand in the limit J — > 0. In 
this limit the evolution operator is site factorizablc and 
given by JJ. e l£/ <f ni (™ i-1 W 2 . This implies periodic time 
dependence, since the operator n,(ni — l)/2 takes integer 
values for any Fock state. An arbitrary initial wave func- 
tion therefore revives entirely after times t„ = 2nn/Uf 
where n is integer (H is set to unity throughout). A 
non vanishing hopping matrix element J greatly com- 
plicates the time evolution, and leads to a relaxation of 
the oscillations. In the following we investigate this sit- 
uation using numerical methods and then interpret the 
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FIG. 1: (Color online) Short time behaviour of the nearest- 
neighbor correlation functions for U% = 2 J and Uf = 40 J. 
(a) t-DMRG and ED results for ID chains. For compari- 
son the thin full line without symbols shows a relaxation for 
Uf — 80 J. (b) ED results for 2D square lattices. In both 
cases there are clear oscillations with a period 2ir/Uf (dotted 
vertical lines for Uf = 40 J), which relax on a time scale of 
1/J. Revival of the correlations at ~ 1.5/ J in (b) is due to 
finite size effects. The two insets show the Fourier transform 
of the oscillations. The main weight lies in a broad band at 
lj m Uf, with some smaller bands at higher multiples of Uf. 



results within a tractable effective model. We implement 
the dynamic transition from the superfluid to the Mott- 
insulating regime at unit filling (n = 1) by a sudden 
quench of the interaction U at fixed hopping J. Thus 
the initial "superfluid" wavefunction is the ground state 
of the Hamiltonian Hi with U — U% < U c , and it evolves 
subject to the Hamiltonian Hj with Uf in the Mott in- 
sulator regime or close to the critical point. The time 
evolution of the many-body wavefunction is computed 
by exact diagonalization (ED) based on Lanczos-type 



methods (171 . |18| and by the adaptive time-dependent 



density matrix renormalization group method (adaptive 
t-DMRG) [HHELi. The ED is used to study ID and 
2D systems with up to 18 sites, while the adaptive t- 
DMRG is used for ID systems with up to 64 sites keeping 
up to 200 DMRG states. For computational reasons the 
Hilbcrt space on each site was truncated at high occupa- 
tion numbers. In ED (t-DMRG) we typically kept up to 
4 (9) bosons per site. 

For a sudden quench deep into the Mott-insulating 
regime we find that a partial revival of the wave-function 
survives. This can for example be seen in the off-diagonal 
correlation functions (b^(x, t)b(0,t)) which display oscil- 
lations with a period set by 1-n/Uf. The oscillations relax 
to a quasi-steady state on a time scale ~ 1/J. We find 
that this relaxation time, as well as the value of the corre- 
lation reached in the quasi-steady state, are independent 
of Uf for sufficiently large Ut. Fig. Q] shows an example 
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FIG. 2: (Color online) Decay of the correlations (&J&i+ r ) with 
distance r after a quench from Ui = 2 J to Uf = 40 J. Squares 
(ED) and triangles (t-DMRG) show the averaged results of 
the quasi steady state. The average value is determined fitting 
a linear function to the results between t\ ~ J -1 and ti = 
20 J -1 . Diamonds show equilibrium QMC results at finite 
temperature (see text for details), and the filled and open 
circles display the T = correlations in the ground state for 
Ui and Uf, respectively. 



of the evolution of nearest neighbor correlations. A sim- 
ilar evolution is observed for longer range correlations. 
Later we use a simple model valid at strong coupling to 
show that the relaxation time is related to the existence 
of a quasi-particle band of width ~ J around an energy 
U f . The short range correlations oscillate with all the 
frequencies in this band, and therefore they dephase af- 
ter a time scale of the order of the inverse band width. 
Indeed a Fourier decomposition of the oscillations (inset 
of Fig. [T|) reveals this band, as well as weaker contribu- 
tions from higher multiples of Uf. The amplitude of the 
oscillations with frequencies of higher multiples depends 
strongly on the particle distribution of the initial state. 

We turn to investigate the nature of the quasi-steady 
state reached by the off-diagonal correlations. The gen- 
eral expectation in non-integrable systems, is that the 
correlations relax to thermal equilibrium. The temper- 
ature at this equilibrium is set by the internal energy 
imposed on the system by the initial conditions. Inter- 
estingly, in spite of the non-integrability of the BHM we 
find two different regimes of behavior depending mostly 
on the magnitude of the interaction strength in the fi- 
nal state. When Uf is very large the correlations in the 
steady state bear strong memory of the initial state. In 
particular their decay with distance is much slower than 
the corresponding thermal correlations and even slower 
than the ground state correlations at the final point. 
This behavior is shown in Fig. [2l The equilibrium fi- 
nite temperature correlations were calculated using quan- 
tum Monte Carlo (QMC) simulations of the BHM (JTJ 
with U = Uf, where the temperature was determined 
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FIG. 3: (Color online) Left panel: Decay of the correlations 
(b\bi+ r ) with distance r after a quench from Ui — J to Uf = 
4J. Squares show the averaged results for times between ti ~ 
J" 1 and ti = 20 J -1 determined as in Fig. [2] Diamonds show 
equilibrium QMC results at finite temperature (T ss 0.8J) 
and U = 4J, and the filled and open circles display the T — 
correlations in the ground state for Ui and Uf, respectively. 
Right panel: Particle number distribution P(n). The labeling 
is the same as in the left panel. 



by matching the on-site particle distribution P(n). This 
yielded T m « 21.5J and T 2D ~ 23.3J. A completely dif- 
ferent regime is realized when the final state is closer to 
the superfluid transition, i.e. Uf < 6 J. In that case the 
correlations at long times do decay with distance faster 
than the ground state correlations and fit reasonably well 
to correlations at thermal equilibrium as calculated using 
QMC. The good fit of the nearest-neighbor correlations 
implies together with the matching of the particle distri- 
bution that the temperature used for the QMC simula- 
tions corresponds to the energy forced into the system 
by the quench. An example of the correlations in this 
regime is displayed in the left panel of Fig. [3] Contrary 
to the regime of large Uf, the oscillations with period 
2ir/Uf are overdamped. We note however that in this 
regime the correlations do not reach a true steady state 
within the time-scale we are able to simulate reliably. 

In contrast to the rich behavior of the off diagonal cor- 
relations, on-site quantities such as the particle distribu- 
tion P(n) reach a much simpler steady state for all U / 
we considered. For large Uf, P(n) almost does not relax 
from the values in the initial state, in agreement with 
mean- field calculation for the case of high fillings (22|. 
For smaller U f some relaxation occurs, see right panel of 
Fig. [31 In all cases we were able to determine a temper- 
ature where the QMC simulations at Uf reproduce the 
long-time behavior of P(n) to good accuracy. 

The essential results of the calculations in different pa- 
rameter regimes are summarized in Fig. [4J The left panel 
presents a non-equilibrium phase diagram in the plane of 
Ui/ J and Uf/J. In one region (large Uf) off-diagonal 
correlations reach a distinctly non thermal steady state. 
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FIG. 4: (Color online) Left panel: Non equilibrium phase 
diagram in the space of initial and final interaction strengths. 
Two regions are found, one where the steady state is distinctly 
non thermal and the other where correlations do appear to 
thermalize within the numerical error bounds. Squares and 
circles mark points in the respective regions where numerical 
results for a one-dimensional system were obtained. The full 
line follows Ui = Uf. Small arrows mark the equilibrium 
critical value for an infinite one-dimensional system. Right 
panel: Difference between the correlation (b\bi+ r ) at steady 
state and its nominal value for the ground state of the final 
Hamiltonian. Data is shown for r = 1 (circles) and r = 2 
(squares) along the cut Ui/J = 2. The crossover from slower 
than ground state decay (5 > 0) to faster than ground state 
is seen at Uf ~ 6J. 



In the other region, the correlations after some time, are 
well described by the thermal equilibrium results. A 
cut along the line Ui/J = 2 shows a clear crossover of 
the steady state correlations from decaying slower than 
ground state correlations at Uf to faster than ground 
state (as expected of thermal correlations) at Uf ss 6 J. 

We now argue that the different steady state regimes 
may be understood qualitatively based on a simple model 
[23j ■ One can think of the excitations of a Mott insulator 
as particles and holes {p\ and h\) that hop on the lattice 
and may be created and annihilated in pairs. From this 
point of view, the initial state (a commensurate super- 
fluid) is a simultaneous condensate of particles and holes. 
Its time dependence is determined by the effective Hamil- 
tonian in the Mott regime. Neglecting quasi-particlc 
interactions, this Hamiltonian is diagonalizcd by a Bo- 
goliubov transformation, so that Hq = J^w a ^kPl^P ak, 
where f3 ak creates a quasi-particle of the Mott insulator. 

Now consider the time evolution of the momentum 
distribution (nk) = (b\j)k)- In terms of the particles 
and holes, a boson is given roughly by the combination 

~ Pit + h-k- Therefore the component (rik(i)) of the 
momentum distribution can be constructed from quasi- 
particles of the Mott insulator carrying momenta k and 
— k, and it oscillates at a frequency t^. This fact bears 
on the dynamics of spatial correlations (bj +r bj), which 
are the Fourier transform of (riu). Short range correla- 
tions receive contributions from all k components, and 
therefore they dephase at a rate comparable to the full 
quasi-particle bandwidth. Long range correlations on the 
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other hand arc dominated by at small k. 

Within this picture thcrmalization occurs due to quasi- 
particle interactions. In particular the quasi-particle pop- 
ulation equilibrates because of quartic processes of the 
form /3q/?k+q/2/3-k+q/2/3o- Note that quartic terms that 
conserve the quasi-particle number cannot by themselves 
induce thermalization because they do not change the 
non-equilibrium quasi-particle population forced on the 
system by the initial conditions. At the level of Fermi's 
golden rule, the process mentioned above may occur only 
while conserving energy. But this is impossible deep 
in the Mott insulator, as long as the gap A is larger 
than half the quasi-particle bandwidth W ~ AzJ. We 
conclude that thcrmalization should be effectively sup- 
pressed in the regime U J, in agreement with the 
numerical results. Of course there are higher order pro- 
cesses that may still induce thermalization, but appar- 
ently this occurs at time scales much larger than the re- 
laxation time 1/J. On the other hand, as the final state 
approaches the critical point A = 0, quasi-particle in- 
teractions become increasingly effective, which leads to 
rapid thcrmalization. The simple picture outlined above 
is strictly valid only for a dilute quasi-particle population, 
that is for initial state close to the transition. However 
from the numerical results it seems to apply more gener- 
ally. 

The results of the present study leave open interesting 
conceptual questions. Is the non thermal state reached a 
true steady state? Is there a longer time scale associated 
with reaching true equilibrium as in the prethermalized 
states discussed in Ref. 24j, or a scale free relaxation 
analogous to aging phenomena? 

We propose that some of these questions may be ad- 
dressed by experiments with ultra cold atoms on optical 
lattices. The simplest observable to consider is the vis- 
ibility of the interference pattern as defined in Q • For 
example, in a one dimensional homogenous tube, follow- 
ing a quench from f/j = 2 J to Uf = 40 J we expect, based 
on t-DMRG calculation, that the visibility will relax to 
a value of approximately 60%, compared to visibility of 
only 20% in the ground state of the system with U = Ur. 

An obstacle for the interpretation of experiments is the 
existence of additional sources of relaxation, most promi- 
nently the confining potential and the presence of many 
parallel tubes. We performed time dependent Gutzwillcr 
calculations to compare the effect of the confining po- 
tential to that of tunneling. The effect of the tunneling 
dominates, if the energy difference between neigbouring 
site Vo(2j — 1) is smaller than the width of the particle- 
hole energy bands ~ 6zJ, i.e. Vo(2j m — l)/(6zJ) < 1. 
Here Vq is the prefactor of the trapping potential and j m 
is the extension of the condensate. Substituting into this 
analysis the experimental parameters of Q it follows that 



the main source of dephasing in that experiment was the 
confining potential. However it is not difficult to reach 
the regime where the tunneling gives the dominant con- 
tribution. 
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